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We show how to use properties of the vectors which are iterated in the transfer-matrix approach 
to Anderson localization, in order to generate the statistical distribution of electronic wavefunction 
amplitudes at arbitrary distances from the origin of L d_1 x oo disordered systems. For d = 1 our 
approach is shown to reproduce exact diagonalization results available in the literature. In d = 2, 
where strips of width L < 64 sites were used, attempted fits of gaussian (log-normal) forms to the 
wavefunction amplitude distributions result in effective localization lengths growing with distance, 
contrary to the prediction from single-parameter scaling theory. We also show that the distributions 
possess a negative skewness S, which is invariant under the usual histogram-collapse rescaling, and 
whose absolute value increases with distance. We find 0.15 < —S < 0.30 for the range of parameters 
used in our study. 
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I. INTRODUCTION 

The localization model introduced by Anderson! in- 
corporates two basic elements, namely the rules of quan- 
tum mechanics applied to a single-electron, tight-binding 
model Hamiltonian, plus quenched disorder (realized, 
e.g., by assigning random self-energies to lattice sites). 
Its original purpose was to show the existence of a 
disorder- induced transition in three-dimensional systems, 
from metallic (diffusive) to insulating (localized) elec- 
tronic behavior, upon increasing randomness. Over 
the years the model has turned out to exhibit a rich 
variety of physical aspects, many of them highlighted 
by the .(single-parameter) scaling theory of localization 
(SPST)BB. Of particular interest here is the fact that, in 
zero magnetic field and in the absence of spin-orbit cou- 
plings, SPST predicts insulating behavior, for any finite 
amount of disorder, in spatial dimensions d = 1 and 2, 
though in the marginal case d — 2 one has borderline phe- 
nomena such as weak localization. Interest in the Ander- 
son transition has been renewed by reports of metallic be- 
havior in dilute two-dimensional electron-hole systemsa. 
While phenomenological, percolation-based theories have 
been able to-jreproduce experimentally-observed trends in 
some detailu, attempts to reconcile basic theoretical as- 
sumptions to experimental evidence have only met lim- 
ited success soiar. For instance, numerical evidence has 
been producedQ against the idea that electron-electron 
interactions (not included in scaling theory) might play 
a role in driving the two-dimensional transitionP. 

Even when one confines oneself to the original An- 
derson picture of non-interacting electrons in three- 
dimensional lattices, where the existence of a transition is 
not questioned, progress towards extracting reliable nu- 
merical estiawrtes-ft£,pitical quantities has been remark- 
ably hardQflla'BliiJ'EJIla. Recently, systematic consider- 
ation of irrelevant variably i-an-d ■ rupnu li pear corrections 
to single-parameter scaling^lla'^OEJM has helped pro- 



duce results with a fairly reasonable claim to consistently 
narrow error bars. 

It is thus of interest to reexamine the basic methods 
which have been used in the past 20 years, in conjunc- 
tion with SPST, to study the Anderson localization prob- 
lem. A step in this direction has been given in Ref. |2(], 
whose authors obtained wave functions via exact diag- 
onalization, for both one-dimensional and finite, L x L, 
two-dimensional systems. By averaging over randomness, 
they obtained probability distributions of wavefunction 
amplitudes on sites at varying distances from an arbitrary 
origin. Such distributions were compared to predictions 
from SPST; though agreement was good in d = 1, the 
two-dimensional results were in contradiction to the idea 
of a single localization length depending only on disor- 
der intensity: instead, a clear logarithmic increase with 
distance, at fixed disorder, was found from their fits for 
that quantity. 

In Ref. finite-size effects were avoided in d — 2 by 
considering disorder strengths such that the correspond- 
ing localization length, as predicted by single-parameter 
theory, is A < 10 (see, e.g., Refs. and using suitably 
large systems with L = 300. On the other hand, nu- 
merous studies of the Anderson transition are set up on 
quasi-one dimensional geometries, for ease of application 
of transfapnHi«drij&-jtEJ^|_^ Green's functions 

methodsQflaMEJJOllalla't0t£l; extrapolation to bulk be- 
havior (d = 2 or 3 as the case may beXis then performed 
with help of finite-size scaling theoryEil. 

Here we consider TM methods, applied both to strictly 
one-dimensional systems and to strips of a square lattice. 
Traditionally the TM approach has been used to calculate 
Lyapunov exponents (directly related to the localization 
length of SPST)LrE3, and quantities obtainable from such 
exponents, e.g., conductancesE§Eil. In Sec. || we recall 
how wavefunction amplitudes may be estimated in the 
TM context, and illustrate our approach in the simple 
d = 1 case by rederiving the corresponding distributions 
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found in Ref. |2(j. In Sec. [I], an analogous treatment 
is developed for strips of a two-dimensional square lat- 
tice, and numerical results are displayed and discussed. 
Conclusions and final remarks are given in Sec. 0. 



II. METHOD AND ONE-DIMENSIONAL 
ILLUSTRATION 

We consider the site-disordered Anderson model, for 
which the tight-binding Hamiltonian is written as 
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(i) 



(id) 



where the site self-energies e< are independent, identically 
distributed random variables obeying a specified distri- 
bution, denotes nearest-neighbor sites on a regular 
lattice, and the energy scale is set by the hopping matrix 
element, V = 1. Disorder intensity is given by the width 
W of the self-energy probability distribution, taken here 
as rectangular (same as in Ref. [2(1) : 



_ /constant - W/2 < e, < +W/2 

-r (£ij — Sq „ii : „_ \^) 



otherwise 



In the TM approachQ, one considers the Hamiltonian 
Eq. (Q) on a quasi-one dimensional L d ~ x x N system, 
N ^> L. Denoting by k = 1, . . . , N the successive cross- 
sections, and i = 1, .. .L d ~ x the respective positions of 
sites within each cross-section, an electronic wave func- 
tion at energy E is given in terms of its local amplitudes, 
{ciik(E) }, and tight-binding orbitals \ik), as: 

$E = ^a ik (E)\ik) . (3) 
With a corresponding change of notation, Eq. reads: 
H = J2 £ ik\ik)(ik\+ > ( 4 ) 

ik (ik,i'k') 

where (ik.i'k') stands for nearest-neighbor pairs. Apply- 
ing Eq. ([|) to Eq. (]|) gives the recursion relation 

Qj,(fc+i) = (E- Ejk) dik - <H,(k-i) ~ ^ Oi>,k , (5) 

i' 

where %' denotes nearest neighbors of % within the same 
cross-section (the ^-dependence is omitted for clarity). 
In matricial form, 
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where, considering, e.g., periodic boundary conditions 
across the d—1 transverse directions!!], 



Pk = 



'E - £i,fc 
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The (2L d ~ 1 x 2i d " x ) matrix T k 



-1 E 

Pk - 
I 



-1 



Erd-\ 







(7) 

is symplec- 



i = 
the 



tic, that is, its eigenvalues occur in pairs {aj, aj 1 } 
I,..., L d ~ x . As explained at length in Refs. 
matrix product Mjy = Y\ k =i T k gives rise to the eigenval- 
ues exp7i . . . exp7 L d-i, where the ji are the Lyapunov 
characteristic exponents (LCE) for the problem, setting 
the asymptotic divergence of the corresponding eigenvec- 
tors Vi. Being a product of symplectic matrices, Mn also 
has this property, therefore the LCE occur in symmet- 
ric pairs {7i,— 7i}. Attention usually concentrates on 
the LCE of smallest modulus, 7i<*-i, whose inverse gives 
the longest decay length (identified with the localization 
length of the Anderson problem) . 

To see the meaning of the eigenvectors v t , recall that 
the physically acceptable (non-diverging) wave function 
is the one associated with the negative LCE of smallest 
modulus, 7£d-i +1 c3; the corresponding tight-binding am- 

plitudes {a ik } will give information on the shape 

of the electronic wave function of interest (because of the 
symplectic character of the TM, one might equally con- 
sider the inverse of the amplitudes associated to 7^-1; 
however, in practice the amount of calculational effort is 
the same either way) . Little use appears to have been 
made of this property in the context of TM studies of 
Anderson localization, except for a calculation of lateral 
transport properties in layered mediaE3. 

It is important to recall that the site amplitudes are 
not the directly relevant quantities in the TM approach; 
instead, in the quasi-one dimensional systems used here 
the wavefunction decay rate must be defined by com- 
paring the Hjgduli of suitable vectors, each with 2L d ~ 1 
componentsBBEJ. 

We now take strictly one-dimensional systems, and il- 
lustrate how the above ideas work. The matrices T k 

are 2x2, and the relevant LCE is 72; the correspond- 

(2) 

ing wavefunction amplitude at site k is a k . Starting 
with an arbitrary pair of states at neighboring sites, say 
(ao, ai) = (1,1), we first iterate Eq. (O) a number iVjn of 
times, taking care to orthonormalize the resulting vectors 
every iV or tho steps (typically, N in — 100; we have used 
iVortho = lj but other authors have used -/V ort ho — Wrf* - 
parently without noticeable deterioration of resultstatd) . 
With such initialization the starting vectors are rotated 
in Hilbert space towards the asymptotic direction of the 
eigenvectors of M/v. Having done this, we rename the 
current site as the origin. Recalling from Eq. (||) that 
the vectors being iterated involve both ip k and ipk+i, one 
sees that the appropriate quantities to keep track of are 
the b ( k 2) = {[a k 2) ] 2 + [4+i] 2 } 1/2 - One might visualize 
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the process as follows: starting from the pair of site am- 
plitudes (ao,ai), one iteratively obtains the pair (01,02) 
and so on, until (after r iterations of T) one gets the pair 
(a r , a r+ i). This latter is legitimately said to be at a dis- 
tance r from the origin, that is, from the original pair of 
sites. 

We then start to accumulate the products of successive 

(2) 

amplitudes b k at each AT or tho steps ( after orthogonal- 
ization, but before normalization). At distance r from 
the new origin, the (relative) wavefunction amplitude is 
given by 



A(r) = - In 



#(0) 



k=l 



(8) 



where the notation of Ref. |2flis used for ease of compar- 
ison, ip'(k) = {[ipk] 2 + [V'/c+if 2 } 1 / 2 , and the fact that we 
always make AT rtho = 1 has been taken into account. In 
order to generate statistics of the A(r) for a set of dis- 
tances {ri < r-i < ■ ■ ■ < r max }, one iterates the TM for 
Ao > r max steps, collecting data at the specified points; 
to collect the next sample, it suffices to keep iterating for 
another Nq steps, with no need to reinitialize the wave 
functions, and so on. After a total of iVjn + N s Nq it- 
erations of the TM, one has N s samples of the A(r) for 
each distance of interest, The corresponding histograms 
H(A, r) for E = 0, disorder strength W = 1.0, r = 1600, 
3200, and 4800, and A" s = 10 5 are shown in Fig. [jj These 
values of the parameters were chosen in order to enable 
comparison with the exact diagonalization results dis- 
played in Figure 1 of Ref. |2(]. Indeed, one finds excel- 
lent visual agreement between the respective data sets. 
Numerical analysis of the first three moments of our his- 
tograms shows that: (i) for all three distances, they are 
well-fitted by gaussians of the form 



H(A,r) 



1/2 



exp 



or I A 



(9) 



with a ~ 2.08, A ~ 104, which are in rather jgopd agree- 
ment with the expected values from SPSTcJ'EZl in the 
limit A s ,r — > 00, respectively a = 2, A = 105.045/W^ 2 ; 
and (ii) the dimensionless skewness S, defined asE3: 



S 




(10) 



for a distribution with mean (x) and dispersion A, has 
the following values: S = 0.0250, 0.0190, 0.0128 respec- 
tively for r = 1600, 3200, 4800. 

Point (ii) is further indication that the wavefunction 
amplitude distribution indeed approaches a log-normal 
shape (zero skewness), but only as r — ► 00. The approx- 
imate dependence on r may be inferred as S ~ r" 
from the 3 data just quoted. 
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Figure 1: Normalized histograms of occurrence of the loga- 
rithmic decay factor, A(r) of Eq. (^), in d = 1, for E — 0, 
W = 1.0, and distances r as shown. N s = 10 samples were 
collected, for each r (see text). 



III. STRIPS OF A TWO-DIMENSIONAL 
LATTICE 



We now extend our approach to strips of a square lat- 
tice. For a strip of width L, the matrices Tk are 2L x 2L, 
and the relevant LCE is 7l+i; the corresponding wave- 
function amplitudes at column k are <& , i = 1, . . . , L. 



Taking into account the normalization of Eq. 
propriate decay factor here is 



A(r) = 



k=l 



In 



fe+l L 

j=k i=i 



1/2 



(L + l) 



the ap- 



(d=2). 



It must be stressed that Eq. (Wv 



is not meant to imply 

an averaging process over site amplitudes 
marked above, these 2L quantities are not the directly 
relevant ones. Instead, they give the modulus of the 
eigenvector associated to the negative LCE of smallest 
absolute valueEa, whose decay is to be followed. 

We have considered strips of even widths 4 < L < 
64 sites and periodic boundary conditions across, and 
taken the disorder intensity W = 10, in order to make 
contact with analogous results in Ref. |2(]. For this value 
of W. SPST predicts the localization length to-be A ~ 
5.45u. SPST, together with finite-size scaling) 2 -!! would 
imply that 



H(A,r,L,W)=f (A, 



V A 



(12) 
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Figure 2: Effective localization lengths on strips of d — 2 
lattice, fitted from first and second moments of distributions 
to the gaussian form, Eq. (^), for E = 0, W = 10.0. Strip 
widths L and distances r as shown. N B = 5 X 10 4 samples for 
L < 40, 2 x 10 4 otherwise. Dashed line is A = 3.4 (logL) 72 
(see text). Scale is logarithmic on horizontal axis. Inset: 
width of distribution, o, as defined in Eq. ([]). Axis scales 
and symbols as in main Figure. 



In order to infer the d — 2 behavior, one must consider 
the regime r, L > 1, r/L < lEil. Although TM methods 
make it easy to explore long distances (r ^> L) along the 
"infinite" direction, this fact is not directly relevant here, 
as the corresponding regime would be one where strictly 
one-dimensional features emerge. 

In the following, we shall restrict ourselves to r/L = 
1/2 and 1; according to Eq. (p^), for each value of r/L one 
should then be able to collapse all distributions against 
r/A. 

Again, we have examined the first three moments of the 
distributions thus generated. Using only the first two, we 
have fitted data to gaussians in the manner of Eq. (jg), 
for which the effective localization lengths A(r, L) are dis- 
played in Fig. ^. It can be seen that, for given L, the A 
are essentially the same both for r = L/2 and L. This 
shows that crude finite-size distortions do not play a role 
for the ranges of r and L used. On the other hand, simi- 
larly to the findings of Ref. ^ and against the prediction 
of SPST, there is no single value of A to fit all distribu- 
tions; instead, it grows with increasing r. However, our 
result differs from that of Ref. |2(], in that the dependence 
of A is clearly not linear in logL. This should not be seen 
as a direct contradiction, as the quantities under study 
are not identical (as was the case in d = 1 ) : though they 
represent the same physical phenomenon of wavefunction 



decay, they do so in rather different geometries. 

As regards the width of distribution a (again taking 
Eq. (j|) as a starting point) our results, displayed in the 
inset of Fig. B e xhibit numerical values not unlike those 
found in Ref. pp|, in the sense of being consistently smaller 
than the SPST prediction a = 2, but with the same 
order of magnitude. One might expect that, for larger 
r, L the decreasing trend observed for 20 < L < 60 
would stabilize close to a ~ 1.3 quoted in Ref. 20, This, 
however, we have no means to ascertain at present. 

We have not been able to fit the full range of data by 
a single power of either L or logL; assuming, e.g., A ~ 
(logL) 2 , the best result from a nonlinear least-squares 
fit gives x ~ 0.72, corresponding to the dashed line in 
Fig. g. It is evident, from the Figure, that the trend for 
large L is towards an even slower variation. 

At this point, one might speculate that A, as given 
by the gaussian fits, could eventually saturate for larger 
L, at a value which might even be close to the SPST 
prediction. However, we shall now show that the gaussian 
fits themselves become increasingly unable to reflect the 
properties of the wavefunction amplitude distributions. 

Indeed, we have found that the skewness of distribu- 
tions is negative, and increases in absolute value as r, 
L grow. As an example, Fig. || shows the raw data 
for L = 48, r = 24, together with the corresponding 
gaussian fit of Eq. (||), obtained using the first and sec- 
ond moments of the distribution. The effect of negative 
skewness is apparent in that the gaussian approximation 
overshoots the data for large A (i.e., predicts a small am- 
plitude to occur more frequently than observed in fact), 
and undershoots for small A (predicts a large amplitude 
to occur less frequently than observed) . 

Skewness data for the ranges of r and L used here 
are displayed in Fig. ||, together with fits of single-power 
forms, — S ~ L x , for the subsets corresponding respec- 
tively to r = L (dashed line, x ~ 0.29) and r = L/2 
(full line, x ~ 0.25). Despite the large amount of scatter, 
the increasing trend against growing r, L is unmistakably 
present. This means, in turn, that for larger and larger 
systems the distributions become ever less amenable to 
fitting by gaussians, as predicted by SPST. This would 
remain true even if a hypothetical saturation should oc- 
cur for values of r and L larger than those investigated 
here. 

Recall, from Eq. (|Io|), that skewness i*-,invariant un- 
der the usual histogram-collapse rescalingEl, H s (A s ,r) = 
H(A, r) nar / \{r) , where the shifted variable is A s = 
[A — r/A(r)] /y/ar/\{r). Therefore, this is a legitimate 
extra parameter to characterize the distributions. Sim- 
ilar results were found experimentally, for the conduc- 
tance distribution in quasi-one dimensional gold wireso. 

Negative skewness of wavefunction amplitude distribu- 
tions works in the same way as (in the limited context 
of gaussian fits) does the finding that A(r) increases with 
r: both contribute to a slower decay of electronic wave 
functions, compared with the constant-A, zero-skewness, 



5 




A 



Figure 3: Normalized histogram of occurrence of A(r) of 
Eq. (0) for E = 0, W = 10.0, L = 48, r = 24 (full line). 
Skewness=— 0.288. Dashed line: Gaussian fit (Eq. (||)), us- 
ing first and second moments of distribution. Vertical scale 
is linear in (a), and logarithmic in (b), the latter in order to 
emphasize discrepancies between data and fit at the extremes. 
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SPST picture. Of course, the evidence just presented is 
not enough to argue that there must be a localization- 
delocalization transition in d — 2; the idea that this is the 
borderline dimensionality, as predicted by SPST, most 
likely holds true. Nonetheless, we have shown robust 
evidence for deviations from SPST in d — 2, whose con- 
sequences still have to be worked out in full, 



IV. CONCLUSIONS 



We have made use of suitable properties of the vec- 
tors which are i^pSit^d in the TM approach to Ander- 
son localizatiorJ3'ut3E3, in order to generate the statisti- 
cal distribution of electronic wavefunction amplitudes at 
sites of x oo disordered systems. We have consid- 

ered d = 1 (for which our approach is shown to repro- 
duce the exact diagonalization results of Ref. ^), and 
d = 2. In the latter case, since the L x oo geometry of 
our systems differs from that (L x L) of Ref. g(| a perfect 
match is not to be expected; however, some basic physical 
properties are found to hold for both cases. In particu- 
lar, attempted fits of gaussian (log-normal) forms to the 
wavefunction amplitude distributions result in effective 
localization lengths growing with distance, contrary to 
the SPST prediction. We have gone further, and shown 
that the distributions possess a negative skewness, which 
is invariant under the usual histogram-collapse rescaling, 
and increases with distance (at least for the range of pa- 
rameters used in our study). 

Such deviations from the expected behavior are evi- 
dence of slower decay of electronic wavefunctions than 
predicted by SPST; it still must be worked out whether 
or not some phenomena specific to d = 2, such as weak lo- 
calization, or the recently-observed metal-insulator tran- 
sition, in dilute two-dimensional dilute electron-hole sys- 
temscl, carry the fingerprints of the anomalies reported 
here. We expect that the present results may motivate 
further work along these lines. 
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